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ABSTRACT 

We carry out three-dimensional, high resolution (up to 1024^ x 256) hydrodynamic simulations of the evolu- 
tion of vortices in vertically unstratified Keplerian disks using the shearing sheet approximation. The transient 
amplification of incompressible, linear amplitude leading waves (which has been proposed as a possible route 
to nonlinear hydrodynamical turbulence in disks) is used as one test of our algorithms; our methods accurately 
capture the predicted amplification, converges at second-order, and is free from aliasing. Waves expected to 
reach nonlinear amplitude at peak amplification become unstable to Kelvin-Helmholtz modes when | W^ax \ ^ ^ 
(where W^ax is the local maximum of vorticity and ^} the angular velocity). We study the evolution of a power- 
law distribution of vorticity consistent with Kolmogorov turbulence; in two-dimensions long-lived vortices 
emerge and decay slowly, similar to previous studies. In three-dimensions, however, vortices are unstable to 
bending modes, leading to rapid decay. Only vortices with a length to width ratio smaller than one survive; in 
three-dimensions the residual kinetic energy and shear stress is at least one order of magnitude smaller than 
in two-dimensions. No evidence for sustained hydrodynamical turbulence and transport is observed in three- 
dimensions. Instead, at late times the residual transport is determined by the amplitude of slowly decaying, 
large-scale vortices (with horizontal extent comparable to the scale height of the disk), with additional contri- 
butions from nearly incompressible inertial waves possible. Evaluating the role that large-scale vortices play in 
astrophysical accretion disks will require understanding the mechanisms that generate and destroy them. 

Subject headings: accretion, accretion disks - solar system: formation 



1. INTRODUCTION 

Many aspects of the structure and evolution of astrophysi- 
cal accretion disks are controlled by the mechanism for an- 
gular momentum transport. When the disk is dynamically 
coupled to a weak magnetic field, the magnetorotational in- 
stability (MRI, e.g., Balbus & Hawley 1998) generates MHD 
turbulence which is a vigorous source of angular momentum 
transport (for a recent review, see Balbus 2003). However, 
in neutral or very weakly ionized disks (e.g., proto-planetary 
disks, and dwarf nova systems in quiescence) the magnetic 
field may be decoupled from the gas, and the MRI may not 
operate (Salmeron & Wardle 2005; Fromang, Terquem, & 
Balbus 2002; Gammie & Menou 1998). It these cases, it is 
important to investigate purely hydrodynamical mechanisms 
for the transport of angular momentum. 

Although Keplerian shear flows are linearly stable accord- 
ing to the Rayleigh criterion, it has long been supposed that 
there exists some mechanism that produces turbulence and 
transport via nonlinear instabilities at high Reynolds num- 
ber. To date, however, no mechanism for such instabilities 
has been explicitly demonstrated, and no evidence for self- 
sustained hydrodynamical turbulence has been found in di- 
rect numerical simulations of Keplerian disks at the same 
Reynolds number needed to capture known nonlinear insta- 
bilities in Cartesian shear flows (Balbus, Hawley, & Stone 
1996; Hawley, Balbus & Winters 1999). Even if nonlinear in- 
stabilities exist in Keplerian shear flows, Lesur & Longaretti 
(2005) argue that current simulations already show the associ- 
ated transport would be negligable. Recently, new interest has 
focused on the role that transient growth of incompressible, 
leading waves might play in generating hydrodynamical tur- 
bulence (Chagelishvili et al. 2003; Umurhan & Regev 2004; 
Yecko 2004). The amplification factor of incompressible (vor- 
tical) waves is of order {kx^o/kyf, where kx^o and ky are the 



initial wavenumbers in the radial and azimuthal directions re- 
spectively (Johnson & Gammie 2005a, hereafter JGa). For 
leading waves that are initially tightly wound | kx^o \ /ky^ 1 , 
large amplification factors are possible, potentially leading to 
nonlinear amplitudes. This has been suggested as one step in 
the route to self- sustained hydrodynamical turbulence in disks 
(Afshordi et al 2005; Mukhopadhyay et al. 2005; but see Bal- 
bus & Hawley 2006). One purpose of this paper is to inves- 
tigate the effect of transient growth of leading waves in fully 
three-dimensional (3D) compressible hydrodynamical disks. 

On the other hand, there is continued interest in the role 
that vortices play in disks. There is particular interest in the 
context of proto-planetary disks, where it has been suggested 
that vortices play a role in trapping dust and aiding in grain- 
growth (e.g.. Barge & Sommeria 1995; Adams & Watkins 
1995; Johansen, Anderson, & Brandenburg 2004). In addi- 
tion, if vortices produce outward angular momentum transport 
they may play an important role in the structure and evolution 
of proto-planetary disks as well as other weakly ionized ac- 
cretion disks (e.g., Li et al. 2001; Barranco & Marcus 2005, 
hereafter BM; Johnson & Gammie 2005b, hereafter JGb). A 
variety of authors have studied the evolution of random vor- 
ticity in both 2D and 3D, and in both incompressible and com- 
pressible gases. Long-lived, anticy clonic (with rotation oppo- 
site to the background shear) vortices have been found in 2D 
Keplerian disks (e.g., Godon & Livio 1999, 2000; Umurhan 
& Regev 2004; JGb). The work of JGb is particularly relevant 
to this study in that they considered the evolution of a spec- 
trum of vorticity consistent with Kolmogorov turbulence in 
fully compressible 2D simulations. In every case, they found 
the kinetic energy and transport decays from the initial condi- 
tions (sustained turbulence was not evident), although the rate 
of decay was affected by factors such as the horizontal ex- 
tent of the domain. These authors argued that compressibility 
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significantly increases the outward angular momentum trans- 
port during the decay. Fully 3D simulations of the evolution 
of vorticity in stratified disks in the anelastic approximation 
have been presented by BM; they found vortices survive (and 
in fact, are produced) in the upper regions of the disk above 
a few scale heights. In the midplane, three-dimensional in- 
stabilities destroy vortices and reduce transport. Here we ex- 
tend the work of JGb and BM by studying vortices in fully 
compressible unstratified disks in 3D. We will report on 3D 
studies of stratified disks elsewhere. 

This paper is organized as follows. In Section 2 we describe 
the numerical methods and shearing sheet model used in our 
simulations. As a test, we study the evolution of both in- 
compressible and compressible planar shearing waves in two- 
dimensions in section 3; we point out that at large amplitude 
vortical waves become Kelvin-Helmholtz unstable. In section 
4 we discuss the evolution of random vortical perturbations in 
both 2D and 3D. We summarize and conclude in section 5. 

2. METHOD 

To study the local hydrodynamics of Keplerian accretion 
disks, the shearing sheet approximation is adopted (Hawley, 
Gammie, & Balbus 1995), that is the equations of motion are 
solved in a Cartesian frame centered at a radius in the disk, 
with linear extent much less than /?o, and which corotates at 
angular velocity Vt{R{)) . In this frame, the Cartesian x coor- 
dinate corresponds to radius and y to azimuth. The hydrody- 
namic equations are expanded to lowest order in H/Rq in this 
frame (H = Cs/^ is the local scale height of the disk, Cs is the 
isothermal sound speed), giving 

^ + V-(pv) = 0, (1) 

1 

at p 

P = c]p , (3) 

where q = -d\x\VL/ dXnRis the shear parameter. For Keplerian 
disks ^ = 3/2. An isothermal equation of state is used through- 
out this paper. In Euler's equation ([2|, the term -2Vt x v is the 
Coriolis force in the rotating frame, and 2qQ^xx represents 
the effect of tidal gravity. The latter is added as the gradient 
of an effective potential, which combined with the fact we use 
a numerical method which solves the equations in conserva- 
tive form (see below), guarantees the net radial and azimuthal 
momentum is conserved. Vertical gravity is neglected for the 
unstratified disks studied here. 

We have omitted an explicit viscosity in our ideal hydrody- 
namics formalism. We show in the next section that the nu- 
merical dissipation inherent in our scheme acts to damp short 
wavelength perturbations as desired (in particular, it prevents 
aliasing of trailing into leading waves), while having little ef- 
fect on well-resolved modes (more than 16 grid points per 
wavelength). It is difficult to define an effective Reynolds 
number for our calculations, since the numerical dissipation 
is a steep function of resolution. We show below the ampli- 
tude errors in the evolution of shearing waves converges at 
better than second-order with our methods. 

Our calculations begin from an equilibrium state consisting 
of constant density p = po, pressure p = pQ = c^pQ, and purely 
azimuthal velocities corresponding to Keplerian rotation v = 
vo = -qftxy. We choose Cg = ft = 10~^ (which implies H =1), 
and po = I. Our computational domain is of size x Ly x 



Lz, with Lx = Ly = 4H and Lz= H. In 2D runs, our domain 
represents the equatorial plane and is of size LxX Ly = AH x 
AH (except in some 2D tests where Lx = Ly = 0.5H, see §§3.2 
and 3.3). Note the difference in the azimuthal velocity across 
the horizontal domain (6c^) exceeds the sound speed. JGb 
found that in the evolution of random vortical perturbations in 
2D, the amplitude and rate of decay of the resulting transport 
depends on the size of this velocity difference (and therefore 
the size of the domain). We adopt the fiducial model studied 
by JGb and expect, as in 2D, that the effect of compressibility 
will be reduced in smaller domains than studied here. 

Our simulations use the 3D version of the Athena code 
(Gardiner & Stone 2005; 2006). Athena is a single-step, di- 
rectionally unsplit Godunov scheme for hydrodynamics and 
MHD which uses the piecewise parabolic method (PPM) 
for spatial reconstruction, a variety of approximate Riemann 
solvers to compute the fluxes, and the corner transport up- 
wind (CTU) unsplit integration algorithm of Colella (1990). 
The calculations presented here use fluxes computed using 
the two-shock approximation for isothermal hydrodynamics. 
The algorithms are second order accurate for all hydrody- 
namic and MHD wave families (Gardiner & Stone 2006); we 
demonstrate better than second-order convergence for plane 
waves in the shearing sheet in §3, indicating that our im- 
plementation of the shearing sheet boundary conditions in 
Athena is at least second-order accurate. 

The computational domain is divided into a grid of Nx x 
Ny X Nz zones. We choose Nx = Ny = N and A^^ = N/A. We 
study the convergence of our numerical results using grids of 
size = 32 to = 2048 in two-dimensions, and = 128 to 
N = 1024 in three-dimensions. 

Important diagnostics of the flows considered here are the 
dimensionless angular momentum flux in fluctuations, de- 
fined as 



{P^xVy)/(pOC,) 



(4) 



where ( ) denotes a volume average, and Vy = Vy + qVtx is the 
difference in the azimuthal velocity compared to Keplerian ro- 
tation. Tests of our numerical methods using epicyclic waves 
show the mean radial flux {pVx) remains zero over hundreds 
of shearing times, thus there is no need to correct equation (4) 
for any mean radial motion. The volume averaged perturbed 
kinetic energy density, including only motions in the horizon- 
tal plane, is 

E^ = {\p{vl + vf)) . (5) 

Finally, the z-component of vorticity 

dvy dvx 



Wz = 



dx dy 



(6) 



is an important diagnostic. In 2D incompressible flows, it is 
exactly conserved (i.e., the Lagrangian derivative dW^/dt is 
zero). In 2D, smooth, compressible flows, the potential vor- 
ticity Wz = (Wz + 2Q) /pis conserved (JGa; the extra 2ft term 
comes from the Coriolis force in the shearing sheet formal- 
ism) if the flow is adiabatic and isentropic. However, in 2D 
flows with shocks both vorticity and potential vorticity can be 
generated through shock-curvature and shock- shock interac- 
tions (Kevlahan 1997). In 3D, can evolve freely. We track 
both Wz and Wz as a diagnostic of our methods, and of the 
effects of shocks. 

3. TWO-DIMENSIONAL EVOLUTION OF SHEARING WAVES 
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Fig. 1. — Evolution of a linear amplitude, compressible, shearing wave, left: The analytic solution (heavy solid line) for the wave amplitude is compared to 
numerical simulations at various grid resolutions N x N, where ^ is the number of grid points per wavelength in the x-direction in the initial conditions, right: 
LI error norm in the j-component of the velocity summed over <Qt <6 versus resolution. The error is converging at N~^-^ . 

We begin by using the linear theory of plane waves in the 
shearing sheet as a test of our numerical algorithms. Our nu- 
merical methods then allow us to investigate what happens 
when leading waves are amplified to nonlinear amplitudes. 
We present the evolution of random vorticity fields in §4. 

In the two-dimensional (x-y) shearing sheet, velocity per- 
turbations representing plane waves have the form 

Vx = ^Vx(t) Qxp[ikx(t)x-\-ikyy] , Vy-\-qQx = Svy(t) Q^xp[ikx(t)x-\-ikyy] 

(7) 

where 

kx(t) = kx^o-\-qnkyt (8) 

and kx^o and > are constant. The subscript denotes quan- 
tities at the initial time ^ = 0. The background shear causes the 
x-component of the wave vector to evolve linearly with time. 
Solutions can be divided into both compressible (V-\y^O) and 
incompressible ( V • v = 0) forms and treated analytically in the 
short- wavelength (i.e., Hky 1), low-frequency (dt <C Cgky) 
limit (e.g. JGa). We consider each separately below. 

3.1. Compressible plane waves 
The amplitude of compressible velocity perturbations 



The oscillatory behavior of the amplitude of a compress- 
ible shearing wave is accurately traced by the numerical re- 
sults. There is no phase error at any resolution, and at ^ = 16 
grid points per wavelength the peak amplitude of the oscilla- 
tions is 0.96 of the analytic value in the first peak, and this 
value increases rapidly with resolution. Below this resolu- 
tion the waves are smoothly damped. As a quantitative mea- 
sure of the errors, figure [J {right) plots the LI error norm 
in the j-component of the velocity summed over the inter- 
val < < 6 as a function of numerical resolution. The 
slope of the line is -2.5, indicating our numerical solution is 
converging at better than second-order. This is direct evidence 
that our implementation of the shearing sheet boundary condi- 
tions, which are required to capture the evolution of shearing 
waves, is no less accurate than our integration algorithm. 

3.2. Incompressible plane waves 

For vortical (incompressible) shearing waves in the short- 
wavelength, low-frequency limit, the amplitudes of the veloc- 
ity perturbation evolve as 



evolve according to (JGa) 



Svx = 6Vx^{)-^ 



5vy = —^5vx , 



(11) 



5vy + {c]k^ + n^)5vy = 



(9) 



y. If the vortical wave starts as leading, i.e.. 



where k,^ = 2(2 - q)ft^ is the square of epicyclic frequency. 
The solution of ^ can be written in terms of parabolic cylin- 
der functions, see JGa for details. Note that equation ^ is 
valid for all wavelengths and frequencies as long as the per- 
turbed potential vorticity is zero (JGa). 

Figure \T\ (left) plots the evolution of a compressible shear- 
ing wave with initial wavenumbers kx^o = -A(2ti/Lx) and ky = 
211 /Ly with Lx = Ly = 4H, and initial amplitude Svx^o/cs = 
4 X 10~^. Also shown is the amplitude of the wave measured 
from numerical simulations at resolutions ofN = 32, 64, 128, 
and 256 grid points per horizontal dimension. An important 
measure of the numerical resolution for plane waves is the 
number of grid points per initial wavelength Xx^o = 27t/ kx^o in 
the x-direction. 



where = kl + 

kxn < 0, it will swing to trailing at a time I /iq^ky). 

During this time, the amplitude of the x- component velocity 
perturbation is temporarily amplified, attaining a maximum 
value Sv^^max = A5v:,fi at t = t^ax, where 



A=l+{k,fi/kyf 



(12) 



Ax 



IttN 

I ^x,0 I 



(10) 



Thus, A/^ = 32, 64, 128, and 256 corresponds to ^ = 8, 16, 32, 
and 64 grid points per initial wavelength respectively. 



is the peak amplification factor (e.g., Chagelishvili et al. 
2003; Umurhan & Regev 2004; JGa). Note that A is also the 
peak amplification factor of the kinetic energy in the wave. 

Figure [2] (left) plots the evolution of an incompressible 
shearing wave with initial wavenumbers ^^^ o = -'^(^n/Lx) and 
ky = 2(27r /Ly), and initial amplitude Svx^o/cs = 10""^. A smaller 
computational domain is used in this test, with Lx = Ly = 0.5 H, 
to minimize the effects of compressibility. To ensure the nu- 
merical representation of the perturbation satisfied the cen- 
tered difference form of V • v = initially, the plane wave solu- 
tion is used to compute the analytic form for the z-component 
of a velocity potential which satisfies v = V x A, and then 
centered differences of Vx = dAz/dy and V3; = -dA^/dx are 
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Fig. 2. — Evolution of a linear amplitude, incompressible, shearing wave, left: The analytic solution (heavy solid line) for the wave amplitude is compared to 
numerical simulations at various grid resolutions N x N, where ^ is the number of grid points per wavelength in the x-direction in the initial conditions, right: 
LI error norm in the x-component of the velocity summed over <Qt <6 versus resolution. The error is converging at N~^-^ for N > 64. 

X 10"^ trailing waves (e.g., fig. 7 in JGb). Figure [3] shows the evolu- 

tion of two runs with parameters fe,o, ^y) = (-I^tt/L^^, "^Tr/Ly), 
Sv^^o/cs = 10"4 and N=64, 128 to Qt = 100. Note that af- 
ter the original peak in the amplification, the amplitude de- 
cays monotonically, without any further transient growth. We 
postulate that the lack of aliasing in our method is related to 
its dissipation properties. At low resolution (below 16 grid 
points per wavelength), Figure[2t shows our method smoothly 
damps waves. Thus, as trailing waves are sheared, our method 
damps them before their wavelength drops to the grid resolu- 
tion and they are aliased into new leading waves. Note this 
does not mean our methods are overly diffusive, because the 
dissipation in Godunov schemes is a strong nonlinear function 
of resolution. Thus, for resolved waves ^ > 16 the dissipation 
is smaller than ZEUS-like methods. In some sense, the dissi- 
pation properties of Godunov methods are ideal for preventing 
aliasing. 

Using fixed resolution per wavelength ^ = 16, we have 
varied the initial ratio kx^o/ky (ky = Air/Ly is fixed) to ob- 
tain different amplification factors A. To keep the waves 
linear at peak amplification, we choose a very small ini- 
tial amplitude 5vx^{)/cs = 10~^. Using grids of Nx = Ny = 
128, 256, 512, 1024, 2048, we have evolved waves with 
I kx^o/ky 1=4, 8, 16, 32, and 64 respectively, giving peak am- 
plification factors of ^= 17, 65, 257, 1025, and 4097. The 
results are plotted in figure H] (left) (dashed lines) along with 
analytical solutions (solid lines). In each case, the numer- 
ical solution tracks the analytic accurately. As long as the 
initial wave is resolved, our method captures large amplifica- 
tion factors as well as it does small. Figure H] (right) plots 
the same data, but scales the amplitude to the peak value, and 
uses r = qQt + kx^o/ky as the time coordinate. The start time 
of each calculation is donated by a filled triangle. In this case, 
every curve lies atop one another, as expected. 

3.3. Kelvin-Helmholtz Instability of Incompressible waves 

Since with sufficient resolution it is possible to accurately 
follow the amplification of incompressible waves by factors 
of 10^"^, it is of interest to investigate what occurs when the 
wave reaches nonlinear amplitude (5vx > c^) at peak ampli- 
fication. Interestingly, we have found that for some initial 
conditions incompressible plane wave solutions become un- 
stable to Kelvin-Helmholtz modes quickly. Figure [5] shows 
the development of the KH instability in shearing wave with 
ky = 47T/Ly, kx^o /ky = -16, slikI initial amplitude Svx^o /cs = 10~^ 



Fig. 3. — Amplification of the amplitude of an incompressible plane wave 
at with kxfi/ky = -4 and resolutions of 64^ and 128^. Aliasing, if present, 
should occur at qflt = Nx/riy -kx^o/ky; note it is completely absent at both 
resolutions. 



used to compute the initial velocities on the grid. Also shown 
is the evolution of the wave amplitude measured from nu- 
merical simulations at resolutions ofA/^ = 32, 64, 128, and 
256 grid points per horizontal dimension, corresponding to 
^ = 8, 16, 32, and 64 grid points per initial wavelength re- 
spectively. 

The expected amplification factor for the wave is A= 17. 
As in the compressible wave case, we find our numerical so- 
lutions accurately track the analytic solution above ^ = 16 
grid points per wavelength. Below this value the wave is 
smoothly damped. Figure^ (right) plots the LI error norm 
in the x-component of the velocity summed over the interval 
0<Qt <6 versus resolution; once again we find uniform con- 
vergence in the errors at a rate of 2.5 for A/^ > 64, that is better 
than second-order. 

An important property of our numerical solutions is the lack 
of any aliasing. As discussed in Umurhan & Regev (2004) and 
JGb, aliasing is dangerous in the sense that it causes power to 
be artificially injected by trailing waves into leading waves, 
which can subsequently by amplified by the shear to repeat 
the loop. The symptom of aliasing is the reoccurrence of tran- 
sient growth at times ftt = (Nx/ny-kx^o/ky)/ q (where ny = 2 
is the y-direction wave number) during the decay phase of 
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Fig. 4. — Evolution of the amplitude of incompressible (vortical) plane waves with different initial | k^fi/ky \ and fixed numerical resolution ^=16 grid 
points per wavelength, left: Curves labelled by = 128,256,512, 1024 and 2048 have | k^^/ky |= 4,8, 16,32, and 64 respectively. The peak amplification 
is ^ = 17, 65, 257, 1025, 4097. Both numerical (dashed line) and analytic (solid line) solutions are shown, right: Same as left, but plotted in terms of 

T = qQt + kx^o/ky. 



on a 5 12^ grid. Four images of Svx are shown at = 0, 0.5, 1 .0 
and 1.5. The destruction of the plane wave is clearly evident. 
We have confirmed that the presence of the instability is in- 
dependent of the numerical resolution. Runs on a 1024^ grid 
show the same behavior, although since the modes are seeded 
by grid noise in the initial conditions, the details change at 
higher resolution. 

Since the KH instability depends only on the amplitude of 
the shear but not the sign, we should expect it to affect both 
leading and trailing waves. By repeating the evolution shown 
in figure[5]but with /ky = +l6,wQ have confirmed that even 
though the wave amplitude is now decaying, it is still KH un- 
stable. In fact, the onset of the KH instability is not related 
to the transient amplification from leading to trailing, but is 
primarily determined by the initial condition. We now give 
some detailed analysis below. 

It is instructive to use the KH instability criterion in planar 
flows to interpret the evolution of plane waves in the shearing 
sheet, even though this criterion ignores the effect of the Cori- 
olis and tidal gravity forces. In a uniform density medium, a 
shear profile with an inflection point is always unstable, with 
a growth rate n ~ kSv. For shearing vortical waves, we might 
expect instability if the KH growth rate exceeds the angular 
velocity, that is 



(13) 



In fact, this crude criterion seems to apply remarkably well 
to our results. For example, simulations with | kx^o \> 1287r 
and Svx^o = 10~^ are unstable, while those with | kx^o \< 327r 
and dvx^o = 10~^ are stable (note ky = liln/Ly) = Stt in these 
tests, and | Svy^o/Svx^o \=\ kx^o/ky |). By "stable" we mean the 
shearing wave solution survives in both the leading phase and 
the trailing phase. 

However, the KH instability of plane shearing waves is not 
related to shear amplification. Note that | ^^^^^3; | in ([T3l) is not 
constant, but evolves with time. The evolution of | ^^^^3; | can 
be described by (if the incompressible plane wave solution is 
stable): 



kxdVy \ = 



Svx^oky I A 



1-- 



1 



(14) 



l + (kx^o/ky-\-q^t) 

Thus I ^^^^^3; I starts from the initial value | Svx^ky \ (^-1), 
decreases in the leading phase and becomes zero at ^max; after 
that, it increases again in the trailing phase and approaches 



the limit of | Svx^oky \ A. Therefore although during the shear 
amplification process, the perturbation amplitude Svx is am- 
plified, the instability condition is not violated. In fact, the 
derivation of criterion ([T3]) only considers azimuthal velocity 
shear. If one further includes radial shear, criterion ([T3]) will 
be modified. The KH instability growth rate is directly related 
to the local maximum of vorticity (e.g., Drazin & Reid 1981), 
i.e., n^ \ Wjnax 1 = 1 kxSvy-kySvx 1 = 1 Svx^oky I A, the latter equal 
sign according to the analytical solution (fill) . Therefore cri- 
terion ([T3]) is modified into 



Wrr. 



SVx^oky 



A>n 



(15) 



Hence the KH instability is determined by the initial condi- 
tions. However, due to the subtlety that the shear is not steady 
and the position of maximum vorticity is changing, the trail- 
ing phase might have a slight preference of KH instability be- 
cause I kx I increases during the trailing phase and the wave 
becomes more and more quasi- steady. 

Figure [6] (/^/O shows the evolution of the volume averaged 
kinetic energy in both leading and trailing waves with param- 
eters Svx^o/cs = 10"^ ky = An/Ly, kx^o/ky = t32 on a 1024^ 
grid along with the expectation from linear theory. For the 
parameters adopted above, peak amplification should occur at 
= 21.3 for the trailing case. For both leading and trailing 
waves, the numerical solution begins to diverge significantly 
from the analytic curve after fit ^2 of evolution. The kinetic 
energy rapidly decays after the KH instability saturates. The 
linear growth rate of the KH instability can be computed from 
the time evolution of the amplitude of the numerical versus 
analytic plane- wave solution; we measure a growth rate of 
about 0.2 I Wmax |, as expected for the fastest growing mode 
(Drazin & Reid 1981). 

Figure[6](ng/z0 shows the evolution of the volume averaged 
angular momentum transport a in both leading and trailing 
waves with the same parameters as in Figure [6^. The KH 
instability causes the value of a to (1) change sign, and (2) 
increase by a factor of about two compared to the expectation 
of linear theory. Note however that at late times a drops to 
near zero, that is there is no evidence for sustained transport. 
Thus, although it might appear attractive to suggest the KH 
instability is a route to turbulence and sustained transport, our 
simulations do not support this idea. We discuss this further 
in §5. 

We suggest that the reason previous studies (Umurhan & 
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Fig. 5. — KH instability in a vortical plane wave. Snapshots of the x-component of the velocity are shown. Parameters are = 512, kxfi/ky = -16, and 
Svx^o/cs = 10~^. Images are at intervals of 0.5Q~^ from the beginning of the run. The planar shearing wave is seriously distorted after only 1 shear time (Qt = 1). 



Regev 2004; JGb) did not find destruction of shearing vortical 
waves is that they used only small initial perturbation ampli- 
tudes and intermediate amplification factors (small | kx^o/ky |), 
which did not satisfy the KH instability condition equation 

4. EVOLUTION OF VORTICES IN TWO- AND THREE-DIMENSIONS 

We now consider the evolution of random vortical (in- 
compressive) velocity perturbations in both two- and three- 
dimensions. Our study closely parallels the work of JGb. We 
draw the initial velocity perturbations from a Gaussian ran- 
dom field with a 2D Kolmogorov power spectrum | S\ p~ 
^-8/3 apply a cut-off to the spectrum at kj^m = 7r/H and 
^max = 32/:min. The amplitude of the random perturbations is 
characterized by a = (| (5v p)^/^, we use a/c^ = 0.8 through- 
out. The effect of varying a was explored in JGb; smaller 
values lead to less angular momentum flux (see below). To 



ensure the perturbations are initially incompressible, we ac- 
tually compute the z-component of a velocity potential 
with the appropriate power spectrum, and then compute the 
velocity perturbations from v = V x A. To study the effect of 
numerical resolution, we have found it critical to compute the 
power spectrum in Fourier space once and for all, and then use 
this same spectrum to generate the initial conditions at every 
resolution. This is only possible because we cut-off the power 
spectrum 3.tk = kj^ax, so there is zero power in the additional 
high-^ modes that can be represented with higher spatial res- 
olution. 

As before, the simulations all use po = 1, Cs = 10~^ , Q = 10~^ 
and a box size of Lx = Ly = AH in 2D (and Lz = H in 3D). 

4.1. Random Vorticity in Two -Dimensions 

We have run simulations with resolutions from A/^ = 128 up 
io N = 2048. Figure |7] shows snapshots in the evolution of 
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Fig. 6. — Evolution of volume-averaged left: kinetic energy, and right: shear stress normalized by the gas pressure in KH unstable leading and trailing 
incompressible plane waves. The parameters are 5v^^q/cs = 10-3, ky = 47T/Ly, k^^o/ky = t32 and = 1024. Light curves are expected analytic solutions and 
heavy curves are numerical results. 



the vorticity fluctuations SWz = Wz + q^ for the 1024^ case. 
As in previous work, we confirm that cyclonic vortices (ro- 
tating in the same sense as background shear) with positive 
values of SW^ are quickly destroyed, while anticyclonic vor- 
tices with negative value of SW^ survive and decay slowly. 
This results in large vortices with negative SW^ embedded in a 
smooth background with a slightly positive SW^. The volume 
average of the vorticity fluctuations reduces to a line integral 
along the boundaries of our domain. Since we use periodic 
boundaries in y, this integral is a measure of the accuracy of 
shearing sheet boundary conditions. We find (SW^) ~ 10~^^, 
i.e. essentially zero, as expected. 

Compressibility does play an important role in the dynam- 
ics; figure[8]is a snapshot of the density and potential vorticity 
Wz at = 20 for the N = 1024 simulation. Strong, interacting 
shocks which propagating primarily in the radial direction are 
clearly evident in the image of the density. The density con- 
trast is more than a factor of three between the lowest and 
highest density regions. The locations of density jumps as- 
sociated with shocks is not immediately evident in the plot 
of Wz- However, horizontal streaks of positive Wz are visible, 
and are associated with shock intersections. As the shocks 
propagate over vortices, they are refracted, producing curva- 
ture in the shock front visible in the density image. At the 
same time, the vortices are impulsively accelerated at each 
shock passage. Since shocks propagate in both radial direc- 
tions, vortices do not gain any net radial motion from shock 
accelerations, but are perturbed back and forth. 

We plot the time evolution of a and £'k at each resolution 
in figure [21 The data points have been boxcar smoothed over 
an interval At = 10Q~^ (a smaller boxcar size ft At = 5 is used 
for data points prior toQt=lO). Compared with fig. 4 in JGb, 
our 128^, 256^ and 512^ runs closely resemble their 256^, 
512^ and 1024^ runs, and show some "convergence" towards 
higher resolution. However, our higher resolution 1024^ and 
2048^ runs do not follow this trend. The time-averaged value 
of a for the last half of the runs (100 < Qt < 200) is of or- 
der ~ 10~^, consistent with the highest resolution runs in JGb. 
Both a and Ek decay by more than an order of magnitude 
during the evolution, and this decay shows no signs of end- 
ing. We find there can be significant differences at late times 
between simulations that use the same initial power-law spec- 
tra, but different initial amplitudes and phases for individual 
Fourier modes. We attribute this to statistical fluctuations in 



the small number of small k modes between different realiza- 
tions of the same power spectrum. The evolution of small k 
modes depends on mode interactions at large amplitudes, and 
the numerical dissipation rate at linear amplitudes. Since the 
numerical dissipation of features resolved by 16 grid points or 
larger is negligible in our methods, once these small k modes 
reach linear amplitude they decay extremely slowly, and can 
dominate the late-time evolution. Thus, the pattern of large- 
scale modes that emerges and decays at late times will de- 
pend on the initial conditions, and the history of mode inter- 
actions that occur during the evolution. Either way, because of 
the small number of small k modes available, we expect sig- 
nificant statistical fluctuations between different runs at late 
times. Our use of boxcar averaging to smooth the data has 
eliminated the high-frequency fluctuations that dominate both 
a and Ek at late times, with typical amplitudes about twice 
the averaged values in the plots. 

The lack of convergence of a and in the 1024^ and 
2048^ runs bears some discussion. Note that at early times, 
both show the same amplitude at ftt = 50, but at a much higher 
level than all the lower resolution runs. Thereafter, the 2048^ 
run decays to the same level as the other runs whereas the 
1024^ run does not. This evolution seems to be related to the 
precise distribution of large-scale vortices that emerge in the 
flow. At Qt = 50, images of the vorticity show four strong, 
well-defined vortices at the higher resolutions, whereas at 
lower resolution only 2-3 weaker and more diffuse vortices 
are evident. There are also a significant number of small, 
weak vortices at higher resolution. Thus, the differing rates 
of diffusion and merging of vortices at different resolutions 
may help to explain the histories. 

4.2. Random Vorticity in Three-Dimensions 

Using numerical resolutions ofA/^=128,256,512 and 1024, 
we have computed the evolution of random vortical perturba- 
tions in 3D. We use the identical initial velocity field, gener- 
ated by the identical Fourier spectrum, as in the 2D case and 
simply project it uniformly over the third dimension. Thus, 
Wz is independent of z in our initial conditions, and therefore 
we study the evolution of vertical vortex columns in 3D. To 
break symmetry in the third dimension, we add small am- 
plitude, random zone-to-zone vertical velocity perturbations 
with maximum amplitude 0.05c^. 

Figure [TOl shows snapshots in the evolution of SWz taken 



8 SHEN, STONE & GARDINER 




Fig. 7. — Evolution of the perturbed z-component of vorticity SWz = Wz + qfl in 2D starting from a random distribution with a power spectrum consistent with 
Kolmogorov at a resolution of 1024^. The box size is 4H x 4H. The initial state is shown in the upper left panel. Frames are taken in lexicographic order at 
Qt = 0, 10, 20, 60. Only anticy clonic vortices (blue) survive at late time. 

from the A/^ = 512 simulation at times of = 0, 10,20 and 
60. Note the pattern of Wz in the x-y plane at r = is similar 
to the first panel of figure Ui Comparing the evolution to the 
2D case, it is clear the vorticity decays much more rapidly. 
Vertical symmetry is maintained until ftt = 20, at which point 
large fluctuations are present as a function of vertical position 
z. By ftt = 60 the initial vortex columns have disintegrated 
into a complex and intertwined network of filaments that show 
no symmetries. The small scale structure introduced by the 
break up of vortex columns likely is the cause of the rapid 
decay. It is well known that columnar vortices are subject to 
elliptical instabilities (e.g., Kerswell 2002), thus the destruc- 
tion of the vertical vortex tubes observed in figure [TO] is not 
surprising. The breakdown of vortex columns has been re- 
ported in numerical simulations (e.g., BM), where the initial 
configurations of vortices are of known analytical forms and 
the growth rate can be measured for isolated vortices. Due to 



the complicated vortex dynamics in our simulation, i.e., back- 
ground shear, vortices interacting and merging, it is difficult 
to quantify the destruction of vortices in terms of elliptical 
instability. Yet the underlying physics should be similar. The 
elliptical instability should exist for all two-dimensional ellip- 
tical streamlines, where 3 -dimensional Kelvin-mode distur- 
bances become unstable when resonating with the underlying 
strain field (e.g., Kerswell 2002). The result of the elliptical 
instability is to break down the vortex column into small-scale 
structures, which is what we find in our simulations. 

Figure [TT] plots the same boxcar smoothed time evolution 
of a and £'k from the 3D simulations, along with the N = 
2048 2D results (heavy solid curves) for reference. The plots 
demonstrate how much more rapid the decay of stress and 
energy is in 3D compared to 2D. The evolution of stress and 
is a rapid exponential decay (from ftt ~ 10-20) followed 
by a slower power-law decay. The exponential decay phase is 
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Fig. 8. — Snapshots of the density (left) and potential vorticity Wz = (Wz + 
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Fig. 9. — Time evolution of the volume-averaged dimensionless shear 
stress a (upper) and 2-dimensional kinetic energy density (bottom) for the 
2D random vorticity runs. 




(right) at Qt = 20 in the 2D random 1024^ run. 

probably associated with the breakup of the vortices and the 
power-law decay phase is probably associated with the decay 
of 3D hydro turbulence. The residual volume averaged shear 
stress is at least one order of magnitude smaller than that of 
the 2D case. 

Perhaps the most striking aspect of the late time evolution 
of a and Ej^ in figure[TT]are the small amplitude (linear) oscil- 
lations. By Fourier analysing the complete (rather than box- 
car averaged) history of a and £^k, we find a strong peak at 
a characteristic frequency of 03ft. Similar, but less strong, 
frequency peaks are also seen at late times in the 2D runs, 
although at a frequency of 1.5Q. This frequency depends 
on boxsize, a 2D box with dimensions Lx = 2H and Ly = 4H 
has a characteristic frequency of the late time oscillations of 
0.751]. Since the oscillations are such small amplitudes, we 
conclude they must be related to linear modes in the shear- 
ing sheet (Balbus 2003). In 2D, only spiral density waves are 
possible, in 3D both acoustic and nearly incompressible iner- 
tial waves are present. The characteristic frequency expected 
depends on the azimuthal wavenumber m, as well as the ratio 
of the vertical to radial wavenumber for inertial waves (e.g., 
Balbus 2003). Since they are nearly incompressible, well- 
resolved (more than 16 grid points per wavelength) inertial 
waves should decay very slowly in our simulations. The ob- 
served frequencies are consistent with low-m linear waves as 
the origin of the oscillations^. 

Figure [12] shows snapshots of Wz at the same time ftt = 5 
at four different resolutions; = 128,256,512 and 1024. No 
convergence in the spatial structure of Wz is seen with resolu- 
tion. Instead, entirely different structures are evident the same 
time at different resolutions. In the lowest resolution case 
N = 128, the pattern of Wz is very smooth, with the remain- 
ing fluctuations fairly symmetric in z. The amplitude of the 
fluctuations increases with resolution, with the vertical struc- 

^ The characteristic frequency of inertial waves is given by ou — mQ = 
^Jk\/(k\ (e.g., Balbus 2003). So the observed low angular frequency 

oscillation (uo r^lir x 0.317) could originate from inertial waves with angular 
frequency u; = (m + y k'^ / (k'^ +k^))Q ~ 1 .88^2, which implies m could be 1. 
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Fig. 10. — Slices of the z component of vorticity Wz in the 512 x 128 3D random vorticity run. Snapshots are taken at VLt = 0, 10, 20, 60. 



ture remaining symmetric up to the N = 5\2 case. However, 
2ii N = 1024, vertical symmetry is broken, and the vorticity 
has disintegrated into small scale filaments similar to those 
observed in the last panel of fi sure [TOl for the = 512 case. 
The lack of convergence in is quite striking. The higher 
resolution simulations are able to capture higher wavenumber 
modes of the elliptical instabilities that lead to the destruction 
of vertical columns. The presence of these modes leads to de- 
struction of the columns at an earlier time at higher resolution. 
This trend accounts for most of the difference in the structure 
between different resolutions in figure [121 The overall trend 
that dynamical instabilities of the vortex columns destroys the 
vertical symmetry and results in more rapid decay of a and £k 
is ubiquitous at every resolution. 

5. SUMMARY AND CONCLUSIONS 

Using the 3D version of the Athena code (Gardiner & Stone 
2005; 2006), we have carried out hydrodynamical simulations 
of the evolution of both planar waves and random vortical per- 
turbations in unstratified Keplerian disks using the shearing 
sheet approximation. Our results can be summarized by the 



following four points. 

(1) Our numerical methods reproduce the evolution of the 
amplitude of both compressible and incompressible (vortical) 
plane waves predicted by linear theory with negligible dis- 
sipation as long as the numerical resolution is 16 grid points 
per wavelength or larger. At lower resolution, the wave ampli- 
tude is smoothly damped. Significantly, there is no evidence 
of aliasing of trailing into leading waves as they are sheared 
down to the grid resolution. We have demonstrated the am- 
plitude error in plane waves in the shearing sheet converges at 
better than second-order with our methods. 

(2) We have shown that incompressible plane waves be- 
come KH unstable and are destroyed when | W^ax |^ ^, the 
condition that the growth rate of the KH instability exceeds 
the angular velocity. Whether the wave becomes KH unstable 
only depends on the initial vorticity amplitude in the initial 
conditions. In this case the planar shearing wave amplitude is 
damped well before it gets amplified by the transient growth 
mechanism. 

(3) In 2D, the evolution of large amplitude random vorticity 
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Fig. 11. — Time evolution of the volume-averaged dimensionless shear 
stress a (upper) and 2-dimensional kinetic energy density (bottom) for the 
3D random vorticity runs. The heavy solid curve is the result of the 2048^ 
2D run for reference. 

perturbations follows the results reported by others (Godon & 
Livio 1999, 2000; Umurhan & Regev 2004; JGb). In partic- 
ular, coherent, large-scale (horizontal extent larger than the 
scale height), anticy clonic vortices survive and decay slowly. 
At late times the dimensionless angular momentum flux is of 
the order 10~^ averaged over 100 <^}t < 200, consistent with 
JGb. This flux is dominated by the residual motions in the 
large-scale anticy clonic vortices. 
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(4) In 3D, the evolution of large amplitude random vortic- 
ity perturbations is similar to the 2D case, except the elliptic 
instabilities destroy vortex columns whose vertical extent ex- 
ceeds their horizontal extent. This greatly increases the rate 
of decay of kinetic energy and stress. The resulting volume 
averaged shear stress a at late-times is at least one order of 
magnitude smaller than that for the 2D case, and is probably 
associated with linear amplitude, low-m inertial waves that re- 
main in the box, and which decay slowly. 

The destruction of planar vortical waves by the KH insta- 
bility has implications for transient amplification as a means 
to drive strong turbulence. Although this might appear at- 
tractive as a feedback mechanism to drive turbulence (Chage- 
HshviH et al. 2003; Yecko 2003; Afshordi, Mukhopadhyay, 
& Narayan 2005), we find instability destroys the wave and 
results in rapid decay of kinetic energy, rather than re-seeding 
new leading waves. We speculate that most of the energy in 
the vortical wave is converted into compressible modes by the 
KH instability, which are not amplified by shear. Moreover, 
Balbus & Hawley (2006) found that these planar vortical wave 
solutions are actually exact to all orders and cannot serve as a 
route to self- sustained turbulence. 

Our results show no evidence for sustained hydrodynami- 
cal turbulence in the shearing sheet. While there is non-zero 
transport, it is associated with large scale vortices that are in- 
troduced in the initial conditions, or which emerge from mode 
interactions during the nonlinear decay phase. At late times 
the Reynolds stress is oscillatory, which may indicate long- 
lived linear amplitude inertial waves also contribute to the 
time-averaged shear. 

To understand the relevance of large-scale vortices to the 
dynamics of astrophysical disks, it will be important to un- 
derstand how they are generated and destroyed. BM have 
shown that such vortices can be generated in stratified disks, 
but given the importance of compressible waves on the decay 
of vortices (JGb), it is important to investigate vortex produc- 
tion and evolution in fully compressible stratified disks. 
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